{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# relax_static calculation style\n",
    "\n",
    "**Lucas M. Hale**, [lucas.hale@nist.gov](mailto:lucas.hale@nist.gov?Subject=ipr-demo), *Materials Science and Engineering Division, NIST*.\n",
    "\n",
    "## Introduction\n",
    "\n",
    "The relax_static calculation style uses static energy/force minimizations to relax the atomic positions and box dimensions of a system to a specified pressure.\n",
    "\n",
    "### Version notes\n",
    "\n",
    "- 2018-07-09: Notebook added.\n",
    "- 2019-07-30: Description updated and small changes due to iprPy version.\n",
    "- 2020-05-22: Version 0.10 update - potentials now loaded from database.\n",
    "- 2020-09-22: Setup and parameter definition streamlined.\n",
    "\n",
    "### Additional dependencies\n",
    "\n",
    "### Disclaimers\n",
    "\n",
    "- [NIST disclaimers](http://www.nist.gov/public_affairs/disclaimer.cfm)\n",
    "- The minimization algorithm will drive the system to a local minimum, which may not be the global minimum.  There is no guarantee that the resulting structure is dynamically stable, and it is possible that the relaxation of certain dimensions may be constrained to move together during the minimization preventing a full relaxation.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Method and Theory\n",
    "\n",
    "This method uses the LAMMPS minimization plus box_relax commands to simultaneously relax both the atomic positions and the system's box dimensions towards a local minimum.  The LAMMPS documentation of the box_relax command notes that the complete minimization algorithm is not well defined which may prevent a complete relaxation during a single run.  To overcome this limitation, the calculation script continuously restarts the minimization until the box dimensions from one run to the next remain within a specified tolerance.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Demonstration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.1. Library imports"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import libraries needed by the calculation. The external libraries used are:\n",
    "\n",
    "- [numpy](http://www.numpy.org/)\n",
    "\n",
    "- [atomman](https://github.com/usnistgov/atomman)\n",
    "\n",
    "- [iprPy](https://github.com/usnistgov/iprPy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Notebook last executed on 2020-09-22 using iprPy version 0.10.2\n"
     ]
    }
   ],
   "source": [
    "# Standard library imports\n",
    "from pathlib import Path\n",
    "import os\n",
    "import shutil\n",
    "import datetime\n",
    "from copy import deepcopy\n",
    "\n",
    "# http://www.numpy.org/\n",
    "import numpy as np  \n",
    "\n",
    "# https://github.com/usnistgov/atomman \n",
    "import atomman as am\n",
    "import atomman.lammps as lmp\n",
    "import atomman.unitconvert as uc\n",
    "\n",
    "# https://github.com/usnistgov/iprPy\n",
    "import iprPy\n",
    "\n",
    "print('Notebook last executed on', datetime.date.today(), 'using iprPy version', iprPy.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.2. Default calculation setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify calculation style\n",
    "calc_style = 'relax_static'\n",
    "\n",
    "# If workingdir is already set, then do nothing (already in correct folder)\n",
    "try:\n",
    "    workingdir = workingdir\n",
    "\n",
    "# Change to workingdir if not already there\n",
    "except:\n",
    "    workingdir = Path('calculationfiles', calc_style)\n",
    "    if not workingdir.is_dir():\n",
    "        workingdir.mkdir(parents=True)\n",
    "    os.chdir(workingdir)\n",
    "    \n",
    "# Initialize connection to library\n",
    "library = iprPy.Library(load=['lammps_potentials'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Assign values for the calculation's run parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.1. Specify system-specific paths\n",
    "\n",
    "- __lammps_command__ is the LAMMPS command to use (required).\n",
    "\n",
    "- __mpi_command__ MPI command for running LAMMPS in parallel. A value of None will run simulations serially."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "lammps_command = 'lmp_serial'\n",
    "mpi_command = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.2. Load interatomic potential\n",
    "\n",
    "- __potential_name__ gives the name of the potential_LAMMPS reference record in the iprPy library to use for the calculation.  \n",
    "\n",
    "- __potential__ is an atomman.lammps.Potential object (required)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "potential_name = '1999--Mishin-Y--Ni--LAMMPS--ipr1'\n",
    "\n",
    "# Retrieve potential and parameter file(s)\n",
    "potential = library.get_lammps_potential(id=potential_name, getfiles=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.3. Load initial unit cell system\n",
    "\n",
    "- __ucell__ is an atomman.System representing a fundamental unit cell of the system (required).  Here, this is generated using the load parameters and symbols."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "avect =  [ 3.500,  0.000,  0.000]\n",
      "bvect =  [ 0.000,  3.500,  0.000]\n",
      "cvect =  [ 0.000,  0.000,  3.500]\n",
      "origin = [ 0.000,  0.000,  0.000]\n",
      "natoms = 4\n",
      "natypes = 1\n",
      "symbols = ('Ni',)\n",
      "pbc = [ True  True  True]\n",
      "per-atom properties = ['atype', 'pos']\n",
      "     id |   atype |  pos[0] |  pos[1] |  pos[2]\n",
      "      0 |       1 |   0.000 |   0.000 |   0.000\n",
      "      1 |       1 |   0.000 |   1.750 |   1.750\n",
      "      2 |       1 |   1.750 |   0.000 |   1.750\n",
      "      3 |       1 |   1.750 |   1.750 |   0.000\n"
     ]
    }
   ],
   "source": [
    "# Create ucell by loading prototype record\n",
    "ucell = am.load('prototype', 'A1--Cu--fcc', symbols='Ni', a=3.5)\n",
    "\n",
    "print(ucell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.4. Modify system\n",
    "\n",
    "- __sizemults__ list of three integers specifying how many times the ucell vectors of $a$, $b$ and $c$ are replicated in creating system.\n",
    "\n",
    "- __system__ is an atomman.System to perform the scan on (required). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# of atoms in system = 108\n"
     ]
    }
   ],
   "source": [
    "sizemults = [3, 3, 3]\n",
    "\n",
    "# Generate system by supersizing ucell\n",
    "system = ucell.supersize(*sizemults)\n",
    "print('# of atoms in system =', system.natoms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.5. Specify calculation-specific run parameters\n",
    "\n",
    "- __pressure_xx__ gives the xx component of the pressure to equilibriate the system to.\n",
    "\n",
    "- __pressure_yy__ gives the yy component of the pressure to equilibriate the system to.\n",
    "\n",
    "- __pressure_zz__ gives the zz component of the pressure to equilibriate the system to.\n",
    "\n",
    "- __pressure_xy__ gives the xy component of the pressure to equilibriate the system to.\n",
    "\n",
    "- __pressure_xz__ gives the xz component of the pressure to equilibriate the system to.\n",
    "\n",
    "- __pressure_yz__ gives the yz component of the pressure to equilibriate the system to.\n",
    "\n",
    "- __displacementkick__ specifies a multiplier for a random shift of atomic positions to apply prior to relaxation.  This is in length units.\n",
    "\n",
    "- __energytolerance__ is the energy tolerance to use during the minimizations. This is unitless.\n",
    "\n",
    "- __forcetolerance__ is the force tolerance to use during the minimizations. This is in energy/length units.\n",
    "\n",
    "- __maxiterations__ is the maximum number of minimization iterations to use.\n",
    "\n",
    "- __maxevaluations__ is the maximum number of minimization evaluations to use.\n",
    "\n",
    "- __maxatommotion__ is the largest distance that an atom is allowed to move during a minimization iteration. This is in length units.\n",
    "\n",
    "- __maxcycles__ is the maximum number of minimization runs (cycles) to perform.\n",
    "\n",
    "- __cycletolerance__ is the relative tolerance to use in identifying if the lattice constants have converged from one cycle to the next. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "pressure_xx = uc.set_in_units(0.0, 'GPa')\n",
    "pressure_yy = uc.set_in_units(0.0, 'GPa')\n",
    "pressure_zz = uc.set_in_units(0.0, 'GPa')\n",
    "pressure_xy = uc.set_in_units(0.0, 'GPa')\n",
    "pressure_xz = uc.set_in_units(0.0, 'GPa')\n",
    "pressure_yz = uc.set_in_units(0.0, 'GPa')\n",
    "displacementkick = uc.set_in_units(0.00001, 'angstrom')\n",
    "energytolerance = 1e-8\n",
    "forcetolerance = uc.set_in_units(0.0, 'eV/angstrom')\n",
    "maxiterations = 10000\n",
    "maxevaluations = 100000\n",
    "maxatommotion = uc.set_in_units(0.01, 'angstrom')\n",
    "maxcycles = 100\n",
    "cycletolerance = 1e-7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Define calculation function(s) and generate template LAMMPS script(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.1. minbox.template"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('minbox.template', 'w') as f:\n",
    "    f.write(\"\"\"# LAMMPS input script that performs an energy minimization and box relaxation\n",
    "\n",
    "box tilt large\n",
    "\n",
    "<atomman_system_pair_info>\n",
    "\n",
    "change_box all triclinic\n",
    "\n",
    "thermo_style custom step lx ly lz xy xz yz pxx pyy pzz pxy pxz pyz pe\n",
    "thermo_modify format float %.13e\n",
    "\n",
    "compute peatom all pe/atom\n",
    "\n",
    "dump dumpit all custom <maxeval> *.dump <dump_keys>\n",
    "dump_modify dumpit format <dump_modify_format>\n",
    "\n",
    "fix boxrelax all box/relax x <p_xx> y <p_yy> z <p_zz> xy <p_xy> xz <p_xz> yz <p_yz>\n",
    "\n",
    "min_modify dmax <dmax>\n",
    "\n",
    "minimize <etol> <ftol> <maxiter> <maxeval>\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.2. relax_static()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def relax_static(lammps_command, system, potential, mpi_command=None,\n",
    "                 p_xx=0.0, p_yy=0.0, p_zz=0.0, p_xy=0.0, p_xz=0.0, p_yz=0.0,\n",
    "                 dispmult=0.0, etol=0.0, ftol=0.0,  maxiter=10000,\n",
    "                 maxeval=100000, dmax=uc.set_in_units(0.01, 'angstrom'),\n",
    "                 maxcycles=100, ctol=1e-10):\n",
    "    \"\"\"\n",
    "    Repeatedly runs the ELASTIC example distributed with LAMMPS until box\n",
    "    dimensions converge within a tolerance.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    lammps_command :str\n",
    "        Command for running LAMMPS.\n",
    "    system : atomman.System\n",
    "        The system to perform the calculation on.\n",
    "    potential : atomman.lammps.Potential\n",
    "        The LAMMPS implemented potential to use.\n",
    "    mpi_command : str, optional\n",
    "        The MPI command for running LAMMPS in parallel.  If not given, LAMMPS\n",
    "        will run serially.\n",
    "    p_xx : float, optional\n",
    "        The value to relax the x tensile pressure component to (default is\n",
    "        0.0).\n",
    "    p_yy : float, optional\n",
    "        The value to relax the y tensile pressure component to (default is\n",
    "        0.0).\n",
    "    p_zz : float, optional\n",
    "        The value to relax the z tensile pressure component to (default is\n",
    "        0.0).\n",
    "    p_xy : float, optional\n",
    "        The value to relax the xy shear pressure component to (default is\n",
    "        0.0).\n",
    "    p_xz : float, optional\n",
    "        The value to relax the xz shear pressure component to (default is\n",
    "        0.0).\n",
    "    p_yz : float, optional\n",
    "        The value to relax the yz shear pressure component to (default is\n",
    "        0.0).\n",
    "    dispmult : float, optional\n",
    "        Multiplier for applying a random displacement to all atomic positions\n",
    "        prior to relaxing. Default value is 0.0.\n",
    "    etol : float, optional\n",
    "        The energy tolerance for the structure minimization. This value is\n",
    "        unitless. (Default is 0.0).\n",
    "    ftol : float, optional\n",
    "        The force tolerance for the structure minimization. This value is in\n",
    "        units of force. (Default is 0.0).\n",
    "    maxiter : int, optional\n",
    "        The maximum number of minimization iterations to use (default is 10000).\n",
    "    maxeval : int, optional\n",
    "        The maximum number of minimization evaluations to use (default is \n",
    "        100000).\n",
    "    dmax : float, optional\n",
    "        The maximum distance in length units that any atom is allowed to relax\n",
    "        in any direction during a single minimization iteration (default is\n",
    "        0.01 Angstroms).\n",
    "    pressure_unit : str, optional\n",
    "        The unit of pressure to calculate the elastic constants in (default is\n",
    "        'GPa').\n",
    "    maxcycles : int, optional\n",
    "        The maximum number of times the minimization algorithm is called.\n",
    "        Default value is 100.\n",
    "    ctol : float, optional\n",
    "        The relative tolerance used to determine if the lattice constants have\n",
    "        converged (default is 1e-10).\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    dict\n",
    "        Dictionary of results consisting of keys:\n",
    "        \n",
    "        - **'relaxed_system'** (*float*) - The relaxed system.\n",
    "        - **'E_coh'** (*float*) - The cohesive energy of the relaxed system.\n",
    "        - **'measured_pxx'** (*float*) - The measured x tensile pressure of the\n",
    "          relaxed system.\n",
    "        - **'measured_pyy'** (*float*) - The measured y tensile pressure of the\n",
    "          relaxed system.\n",
    "        - **'measured_pzz'** (*float*) - The measured z tensile pressure of the\n",
    "          relaxed system.\n",
    "        - **'measured_pxy'** (*float*) - The measured xy shear pressure of the\n",
    "          relaxed system.\n",
    "        - **'measured_pxz'** (*float*) - The measured xz shear pressure of the\n",
    "          relaxed system.\n",
    "        - **'measured_pyz'** (*float*) - The measured yz shear pressure of the\n",
    "          relaxed system.\n",
    "    \"\"\"\n",
    "    # Build filedict if function was called from iprPy\n",
    "    try:\n",
    "        assert __name__ == pkg_name\n",
    "        calc = iprPy.load_calculation(calculation_style)\n",
    "        filedict = calc.filedict\n",
    "    except:\n",
    "        filedict = {}\n",
    "    \n",
    "    # Get lammps units\n",
    "    lammps_units = lmp.style.unit(potential.units)\n",
    "    \n",
    "    # Get lammps version date\n",
    "    lammps_date = lmp.checkversion(lammps_command)['date']\n",
    "    \n",
    "    # Save initial configuration as a dump file\n",
    "    system.dump('atom_dump', f='initial.dump')\n",
    "    \n",
    "    # Apply small random distortions to atoms\n",
    "    system.atoms.pos += dispmult * np.random.rand(*system.atoms.pos.shape) - dispmult / 2\n",
    "    \n",
    "    # Initialize parameters\n",
    "    old_vects = system.box.vects\n",
    "    converged = False\n",
    "    \n",
    "    # Run minimizations up to maxcycles times\n",
    "    for cycle in range(maxcycles):\n",
    "        \n",
    "        # Define lammps variables\n",
    "        lammps_variables = {}\n",
    "        system_info = system.dump('atom_data', f='init.dat',\n",
    "                                  potential=potential,\n",
    "                                  return_pair_info=True)\n",
    "        lammps_variables['atomman_system_pair_info'] = system_info\n",
    "        \n",
    "        lammps_variables['p_xx'] = uc.get_in_units(p_xx, lammps_units['pressure'])\n",
    "        lammps_variables['p_yy'] = uc.get_in_units(p_yy, lammps_units['pressure'])\n",
    "        lammps_variables['p_zz'] = uc.get_in_units(p_zz, lammps_units['pressure'])\n",
    "        lammps_variables['p_xy'] = uc.get_in_units(p_xy, lammps_units['pressure'])\n",
    "        lammps_variables['p_xz'] = uc.get_in_units(p_xz, lammps_units['pressure'])\n",
    "        lammps_variables['p_yz'] = uc.get_in_units(p_yz, lammps_units['pressure'])\n",
    "        lammps_variables['etol'] = etol\n",
    "        lammps_variables['ftol'] = uc.get_in_units(ftol, lammps_units['force'])\n",
    "        lammps_variables['maxiter'] = maxiter\n",
    "        lammps_variables['maxeval'] = maxeval\n",
    "        lammps_variables['dmax'] = uc.get_in_units(dmax, lammps_units['length'])\n",
    "        \n",
    "        # Set dump_keys based on atom_style\n",
    "        if potential.atom_style in ['charge']:\n",
    "            lammps_variables['dump_keys'] = 'id type q x y z c_peatom'\n",
    "        else:\n",
    "            lammps_variables['dump_keys'] = 'id type x y z c_peatom'\n",
    "        \n",
    "        # Set dump_modify_format based on lammps_date\n",
    "        if lammps_date < datetime.date(2016, 8, 3):\n",
    "            if potential.atom_style in ['charge']:\n",
    "                lammps_variables['dump_modify_format'] = '\"%d %d %.13e %.13e %.13e %.13e %.13e\"'\n",
    "            else:\n",
    "                lammps_variables['dump_modify_format'] = '\"%d %d %.13e %.13e %.13e %.13e\"'\n",
    "        else:\n",
    "            lammps_variables['dump_modify_format'] = 'float %.13e'\n",
    "        \n",
    "        # Write lammps input script\n",
    "        template_file = 'minbox.template'\n",
    "        lammps_script = 'minbox.in'\n",
    "        template = iprPy.tools.read_calc_file(template_file, filedict)\n",
    "        with open(lammps_script, 'w') as f:\n",
    "            f.write(iprPy.tools.filltemplate(template, lammps_variables, '<', '>'))\n",
    "        \n",
    "        # Run LAMMPS and extract thermo data\n",
    "        logfile = 'log-' + str(cycle) + '.lammps'\n",
    "        output = lmp.run(lammps_command, lammps_script, mpi_command, logfile=logfile)\n",
    "        thermo = output.simulations[0]['thermo']\n",
    "        \n",
    "        # Clean up dump files\n",
    "        Path('0.dump').unlink()\n",
    "        last_dump_file = str(thermo.Step.values[-1]) + '.dump'\n",
    "        renamed_dump_file = 'relax_static-' + str(cycle) + '.dump'\n",
    "        shutil.move(last_dump_file, renamed_dump_file)\n",
    "        \n",
    "        # Load relaxed system\n",
    "        system = am.load('atom_dump', renamed_dump_file, symbols=system.symbols)\n",
    "        \n",
    "        # Test if box dimensions have converged\n",
    "        if np.allclose(old_vects, system.box.vects, rtol=ctol, atol=0):\n",
    "            converged = True\n",
    "            break\n",
    "        else:\n",
    "            old_vects = system.box.vects\n",
    "    \n",
    "    # Check for convergence\n",
    "    if converged is False:\n",
    "        raise RuntimeError('Failed to converge after ' + str(maxcycles) + ' cycles')\n",
    "    \n",
    "    # Zero out near-zero tilt factors\n",
    "    lx = system.box.lx\n",
    "    ly = system.box.ly\n",
    "    lz = system.box.lz\n",
    "    xy = system.box.xy\n",
    "    xz = system.box.xz\n",
    "    yz = system.box.yz\n",
    "    if np.isclose(xy/ly, 0.0, rtol=0.0, atol=1e-10):\n",
    "        xy = 0.0\n",
    "    if np.isclose(xz/lz, 0.0, rtol=0.0, atol=1e-10):\n",
    "        xz = 0.0\n",
    "    if np.isclose(yz/lz, 0.0, rtol=0.0, atol=1e-10):\n",
    "        yz = 0.0\n",
    "    system.box.set(lx=lx, ly=ly, lz=lz, xy=xy, xz=xz, yz=yz)\n",
    "    system.wrap()\n",
    "    \n",
    "    # Build results_dict\n",
    "    results_dict = {}\n",
    "    results_dict['dumpfile_initial'] = 'initial.dump'\n",
    "    results_dict['symbols_initial'] = system.symbols\n",
    "    results_dict['dumpfile_final'] = renamed_dump_file\n",
    "    results_dict['symbols_final'] = system.symbols\n",
    "    results_dict['E_coh'] = uc.set_in_units(thermo.PotEng.values[-1] / system.natoms,\n",
    "                                       lammps_units['energy'])\n",
    "                                       \n",
    "    results_dict['lx'] = uc.set_in_units(lx, lammps_units['length'])\n",
    "    results_dict['ly'] = uc.set_in_units(ly, lammps_units['length'])\n",
    "    results_dict['lz'] = uc.set_in_units(lz, lammps_units['length'])\n",
    "    results_dict['xy'] = uc.set_in_units(xy, lammps_units['length'])\n",
    "    results_dict['xz'] = uc.set_in_units(xz, lammps_units['length'])\n",
    "    results_dict['yz'] = uc.set_in_units(yz, lammps_units['length'])\n",
    "    \n",
    "    results_dict['measured_pxx'] = uc.set_in_units(thermo.Pxx.values[-1],\n",
    "                                                   lammps_units['pressure'])\n",
    "    results_dict['measured_pyy'] = uc.set_in_units(thermo.Pyy.values[-1],\n",
    "                                                   lammps_units['pressure'])\n",
    "    results_dict['measured_pzz'] = uc.set_in_units(thermo.Pzz.values[-1],\n",
    "                                                   lammps_units['pressure'])\n",
    "    results_dict['measured_pxy'] = uc.set_in_units(thermo.Pxy.values[-1],\n",
    "                                                   lammps_units['pressure'])\n",
    "    results_dict['measured_pxz'] = uc.set_in_units(thermo.Pxz.values[-1],\n",
    "                                                   lammps_units['pressure'])\n",
    "    results_dict['measured_pyz'] = uc.set_in_units(thermo.Pyz.values[-1],\n",
    "                                                   lammps_units['pressure'])\n",
    "    \n",
    "    return results_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Run calculation function(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_dict = relax_static(lammps_command, system, potential,\n",
    "                            mpi_command = mpi_command,\n",
    "                            p_xx = pressure_xx, \n",
    "                            p_yy = pressure_yy, \n",
    "                            p_zz = pressure_zz,\n",
    "                            p_xy = pressure_xy, \n",
    "                            p_xz = pressure_xz, \n",
    "                            p_yz = pressure_yz,                            \n",
    "                            dispmult = displacementkick,\n",
    "                            etol = energytolerance,\n",
    "                            ftol = forcetolerance,\n",
    "                            maxiter = maxiterations,\n",
    "                            maxeval = maxevaluations,\n",
    "                            dmax = maxatommotion,\n",
    "                            maxcycles = maxcycles,\n",
    "                            ctol = cycletolerance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['dumpfile_initial', 'symbols_initial', 'dumpfile_final', 'symbols_final', 'E_coh', 'lx', 'ly', 'lz', 'xy', 'xz', 'yz', 'measured_pxx', 'measured_pyy', 'measured_pzz', 'measured_pxy', 'measured_pxz', 'measured_pyz'])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results_dict.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Report results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.1. Define units for outputting values\n",
    "\n",
    "- __length_unit__ is the unit of length to display values in.\n",
    "- __energy_unit__ is the unit of energy to display values in.\n",
    "- __pressure_unit__ is the unit of pressure to display values in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "length_unit = 'angstrom'\n",
    "energy_unit = 'eV'\n",
    "pressure_unit = 'GPa'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.2. Print Ecoh and lattice constants of relaxed ucell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ecoh = -4.449999998325555 eV\n",
      "a = 3.51999946361013 angstrom\n",
      "b = 3.5199994969059634 angstrom\n",
      "c = 3.519999505001706 angstrom\n",
      "alpha = 90.0\n",
      "beta =  90.0\n",
      "gamma = 90.0\n"
     ]
    }
   ],
   "source": [
    "print('Ecoh =', uc.get_in_units(results_dict['E_coh'], energy_unit), energy_unit)\n",
    "\n",
    "box = am.Box(lx=results_dict['lx'], ly=results_dict['ly'], lz=results_dict['lz'],\n",
    "             xy=results_dict['xy'], xz=results_dict['xz'], yz=results_dict['yz'])\n",
    "\n",
    "print('a =', uc.get_in_units(box.a / sizemults[0], length_unit), length_unit)\n",
    "print('b =', uc.get_in_units(box.b / sizemults[1], length_unit), length_unit) \n",
    "print('c =', uc.get_in_units(box.c / sizemults[2], length_unit), length_unit) \n",
    "print('alpha =', box.alpha)\n",
    "print('beta = ', box.beta)\n",
    "print('gamma =', box.gamma)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.3. Check final system pressures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pxx = -7.1585342727316e-06 GPa\n",
      "Pyy = -8.1046080506064e-06 GPa\n",
      "Pzz = -8.3345122489377e-06 GPa\n",
      "Pyz = 3.6844223558997003e-10 GPa\n",
      "Pxz = 1.6824819411886998e-10 GPa\n",
      "Pxy = 3.9876281320881e-11 GPa\n"
     ]
    }
   ],
   "source": [
    "print('Pxx =', uc.get_in_units(results_dict['measured_pxx'], pressure_unit), pressure_unit)\n",
    "print('Pyy =', uc.get_in_units(results_dict['measured_pyy'], pressure_unit), pressure_unit)\n",
    "print('Pzz =', uc.get_in_units(results_dict['measured_pzz'], pressure_unit), pressure_unit)\n",
    "print('Pyz =', uc.get_in_units(results_dict['measured_pyz'], pressure_unit), pressure_unit)\n",
    "print('Pxz =', uc.get_in_units(results_dict['measured_pxz'], pressure_unit), pressure_unit)\n",
    "print('Pxy =', uc.get_in_units(results_dict['measured_pxy'], pressure_unit), pressure_unit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.4. Show relaxed atomic configuration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "avect =  [10.560,  0.000,  0.000]\n",
      "bvect =  [ 0.000, 10.560,  0.000]\n",
      "cvect =  [ 0.000,  0.000, 10.560]\n",
      "origin = [-0.030, -0.030, -0.030]\n",
      "natoms = 108\n",
      "natypes = 1\n",
      "symbols = ('Ni',)\n",
      "pbc = [ True  True  True]\n",
      "per-atom properties = ['atype', 'pos', 'atom_id', 'c_peatom']\n",
      "     id |   atype |  pos[0] |  pos[1] |  pos[2]\n",
      "      0 |       1 |  10.530 |  -0.030 |  10.530\n",
      "      1 |       1 |  10.530 |   1.730 |   1.730\n",
      "      2 |       1 |   1.730 |  -0.030 |   1.730\n",
      "      3 |       1 |   1.730 |   1.730 |  10.530\n",
      "      4 |       1 |   3.490 |  10.530 |  10.530\n",
      "      5 |       1 |   3.490 |   1.730 |   1.730\n",
      "      6 |       1 |   5.250 |  -0.030 |   1.730\n",
      "      7 |       1 |   5.250 |   1.730 |  -0.030\n",
      "      8 |       1 |   7.010 |  -0.030 |  -0.030\n",
      "      9 |       1 |   7.010 |   1.730 |   1.730\n",
      "     10 |       1 |   8.770 |  10.530 |   1.730\n",
      "     11 |       1 |   8.770 |   1.730 |  -0.030\n",
      "     12 |       1 |  -0.030 |   3.490 |  -0.030\n",
      "     13 |       1 |  10.530 |   5.250 |   1.730\n",
      "     14 |       1 |   1.730 |   3.490 |   1.730\n",
      "     15 |       1 |   1.730 |   5.250 |  10.530\n",
      "     16 |       1 |   3.490 |   3.490 |  -0.030\n",
      "     17 |       1 |   3.490 |   5.250 |   1.730\n",
      "     18 |       1 |   5.250 |   3.490 |   1.730\n",
      "     19 |       1 |   5.250 |   5.250 |  10.530\n",
      "     20 |       1 |   7.010 |   3.490 |  -0.030\n",
      "     21 |       1 |   7.010 |   5.250 |   1.730\n",
      "     22 |       1 |   8.770 |   3.490 |   1.730\n",
      "     23 |       1 |   8.770 |   5.250 |  -0.030\n",
      "     24 |       1 |  -0.030 |   7.010 |  10.530\n",
      "     25 |       1 |  10.530 |   8.770 |   1.730\n",
      "     26 |       1 |   1.730 |   7.010 |   1.730\n",
      "     27 |       1 |   1.730 |   8.770 |  -0.030\n",
      "     28 |       1 |   3.490 |   7.010 |  -0.030\n",
      "     29 |       1 |   3.490 |   8.770 |   1.730\n",
      "     30 |       1 |   5.250 |   7.010 |   1.730\n",
      "     31 |       1 |   5.250 |   8.770 |  -0.030\n",
      "     32 |       1 |   7.010 |   7.010 |  10.530\n",
      "     33 |       1 |   7.010 |   8.770 |   1.730\n",
      "     34 |       1 |   8.770 |   7.010 |   1.730\n",
      "     35 |       1 |   8.770 |   8.770 |  -0.030\n",
      "     36 |       1 |  10.530 |  10.530 |   3.490\n",
      "     37 |       1 |  -0.030 |   1.730 |   5.250\n",
      "     38 |       1 |   1.730 |  -0.030 |   5.250\n",
      "     39 |       1 |   1.730 |   1.730 |   3.490\n",
      "     40 |       1 |   3.490 |  -0.030 |   3.490\n",
      "     41 |       1 |   3.490 |   1.730 |   5.250\n",
      "     42 |       1 |   5.250 |  -0.030 |   5.250\n",
      "     43 |       1 |   5.250 |   1.730 |   3.490\n",
      "     44 |       1 |   7.010 |  -0.030 |   3.490\n",
      "     45 |       1 |   7.010 |   1.730 |   5.250\n",
      "     46 |       1 |   8.770 |  -0.030 |   5.250\n",
      "     47 |       1 |   8.770 |   1.730 |   3.490\n",
      "     48 |       1 |  10.530 |   3.490 |   3.490\n",
      "     49 |       1 |  -0.030 |   5.250 |   5.250\n",
      "     50 |       1 |   1.730 |   3.490 |   5.250\n",
      "     51 |       1 |   1.730 |   5.250 |   3.490\n",
      "     52 |       1 |   3.490 |   3.490 |   3.490\n",
      "     53 |       1 |   3.490 |   5.250 |   5.250\n",
      "     54 |       1 |   5.250 |   3.490 |   5.250\n",
      "     55 |       1 |   5.250 |   5.250 |   3.490\n",
      "     56 |       1 |   7.010 |   3.490 |   3.490\n",
      "     57 |       1 |   7.010 |   5.250 |   5.250\n",
      "     58 |       1 |   8.770 |   3.490 |   5.250\n",
      "     59 |       1 |   8.770 |   5.250 |   3.490\n",
      "     60 |       1 |  -0.030 |   7.010 |   3.490\n",
      "     61 |       1 |  -0.030 |   8.770 |   5.250\n",
      "     62 |       1 |   1.730 |   7.010 |   5.250\n",
      "     63 |       1 |   1.730 |   8.770 |   3.490\n",
      "     64 |       1 |   3.490 |   7.010 |   3.490\n",
      "     65 |       1 |   3.490 |   8.770 |   5.250\n",
      "     66 |       1 |   5.250 |   7.010 |   5.250\n",
      "     67 |       1 |   5.250 |   8.770 |   3.490\n",
      "     68 |       1 |   7.010 |   7.010 |   3.490\n",
      "     69 |       1 |   7.010 |   8.770 |   5.250\n",
      "     70 |       1 |   8.770 |   7.010 |   5.250\n",
      "     71 |       1 |   8.770 |   8.770 |   3.490\n",
      "     72 |       1 |  -0.030 |  10.530 |   7.010\n",
      "     73 |       1 |  10.530 |   1.730 |   8.770\n",
      "     74 |       1 |   1.730 |  -0.030 |   8.770\n",
      "     75 |       1 |   1.730 |   1.730 |   7.010\n",
      "     76 |       1 |   3.490 |  10.530 |   7.010\n",
      "     77 |       1 |   3.490 |   1.730 |   8.770\n",
      "     78 |       1 |   5.250 |  -0.030 |   8.770\n",
      "     79 |       1 |   5.250 |   1.730 |   7.010\n",
      "     80 |       1 |   7.010 |  -0.030 |   7.010\n",
      "     81 |       1 |   7.010 |   1.730 |   8.770\n",
      "     82 |       1 |   8.770 |  -0.030 |   8.770\n",
      "     83 |       1 |   8.770 |   1.730 |   7.010\n",
      "     84 |       1 |  -0.030 |   3.490 |   7.010\n",
      "     85 |       1 |  -0.030 |   5.250 |   8.770\n",
      "     86 |       1 |   1.730 |   3.490 |   8.770\n",
      "     87 |       1 |   1.730 |   5.250 |   7.010\n",
      "     88 |       1 |   3.490 |   3.490 |   7.010\n",
      "     89 |       1 |   3.490 |   5.250 |   8.770\n",
      "     90 |       1 |   5.250 |   3.490 |   8.770\n",
      "     91 |       1 |   5.250 |   5.250 |   7.010\n",
      "     92 |       1 |   7.010 |   3.490 |   7.010\n",
      "     93 |       1 |   7.010 |   5.250 |   8.770\n",
      "     94 |       1 |   8.770 |   3.490 |   8.770\n",
      "     95 |       1 |   8.770 |   5.250 |   7.010\n",
      "     96 |       1 |  -0.030 |   7.010 |   7.010\n",
      "     97 |       1 |  -0.030 |   8.770 |   8.770\n",
      "     98 |       1 |   1.730 |   7.010 |   8.770\n",
      "     99 |       1 |   1.730 |   8.770 |   7.010\n",
      "    100 |       1 |   3.490 |   7.010 |   7.010\n",
      "    101 |       1 |   3.490 |   8.770 |   8.770\n",
      "    102 |       1 |   5.250 |   7.010 |   8.770\n",
      "    103 |       1 |   5.250 |   8.770 |   7.010\n",
      "    104 |       1 |   7.010 |   7.010 |   7.010\n",
      "    105 |       1 |   7.010 |   8.770 |   8.770\n",
      "    106 |       1 |   8.770 |   7.010 |   8.770\n",
      "    107 |       1 |   8.770 |   8.770 |   7.010\n"
     ]
    }
   ],
   "source": [
    "finalsystem = am.load('atom_dump', results_dict['dumpfile_final'],\n",
    "                      symbols=results_dict['symbols_final'])\n",
    "print(finalsystem)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
